The steered molecular dynamics (SMD) simulations mimic force spectroscopy measurement in single molecule manipulation experiment like AFM, but proceed on much shorter time scale.^[26] In a SMD simulation, one terminal end of the protein is fixed, while the other end is coupled to one end of a dummy spring with elastic constant k. When the end of the dummy spring moves away from the protein with constant speed v, its position ξ is a linear function of time

During the pulling process, the force on the end of the protein, f, is calculated from the extension of the dummy spring,

Fig. 2.

Figure Option ViewDownloadNew Window

Fig. 2. Sketch of SMD simulation and the umbrella sampling method. (a) In SMD simulation, one terminus of the protein of interest is fixed, while the other terminus is coupled with a dummy spring with linear elastic constant k. The other end of the spring moves away from the protein with constant velocity v to stretch the protein. (b) In umbrella sampling, in conformational space projected to the extension as the reaction coordinate, the folded state and the unfolded state are located at each end of the coordinate. Between the folded state and the unfolded state, a spring is used to restrain conformations in a small window of conformational space. The system is well sampled in each window, and overlapping between neighboring windows helps to build PMF of the protein.

One of the main challenges in computational biology is the calculation of free energy difference. Here we use umbrella sampling to construct PMF of proteins.^[28] In umbrella sampling, a system is driven to conformations with high free energy by a bias potential along a reaction coordinate.^[29] The conformations are covered by a series of overlapping windows (Fig. 2(b)), in each of which a bias potential, w_i(x), serves to confine the conformational space to achieve a more efficient sampling.^[30] In principle, the bias potentials can have any functional form, but harmonic potentials are generally chosen for their simplicity,

An accurate calculation of PMF requires an adequate number of windows (separate simulations) with enough overlap in conformational space.^[28,31] Based on SMD simulation, a series of conformations were extracted from the SMD trajectory at regular interval of extension. Then each extracted conformation was used as the initial conformation to perform the umbrella sampling simulation in each window. When all simulations were done, the weighted histogram analysis method (WHAM), an extension of Ferrenberg and Swendsens multiple histogram technique, was applied to properly merge the data from all simulations to obtain the global PMF profile.^{[32, 33]}

Here we list detail parameters of SMD and umbrella simulations. Proteins GB1 and gpW with explicit waters were simulated in the NPT ensemble using GROMACS package 5.1.4, with the GROMOS96 53A6 all-atom force field at the temperature of 310 K. A time step of 2 fs was selected. Periodic boundary conditions were applied in all directions. Particle-mesh-Ewald summation was selected for the long-range Coulomb interactions. The cut-off distances for electrostatics and vdW interactions were set at 1.4 nm. The simulations started from an energy-minimized average NMR structure, deposited as entry 1PGA for protein GB1 and 2L6Q for protein gpW in the RCSB protein data bank. The structures of both proteins were solvated with SPC model water before simulations. The pulling velocity of the restraint point in the SMD simulation was 0.001 nm/ps, and the spring constant was 100 kJ/mol/nm². The umbrella simulations were performed in 95 extension windows for GB1 and 60 windows for gpW with spring constant of 500 kJ/mol/nm². Each simulation lasted for 15 ns, and the trajectory from 5 ns to 15 ns was used for the analysis while the first 5 ns simulations were applied for reaching equilibrium.

SMD simulation of GB1 and gpW began with an equilibrated folded structure and stopped when reaching an extension of ∼ 7.1 nm after ∼ 5 ns of simulation. GB1 and gpW show different capabilities to sustain the stretching force in SMD simulation with the same pulling velocity and spring constant.

Figure 3(a) shows the force–extension profiles of the pulling processes of protein GB1. Initially, the force grows linearly while the extension of the protein grows within a limited scale of ∼ 2.69 nm to ∼ 3.20 nm, which indicates that the protein maintains a mechanical stable structure due to the strong force-bearing hydrogen bonds between the first and last β strands. At the force peak with extension of ∼ 3.19 nm, the force drops suddenly and the extension increases to more than 7 nm due to broken of these hydrogen bonds. The corresponding conformation at the force peak is supposed to be the transition state of force-induced unfolding of GB1.

Fig. 3.

Figure Option ViewDownloadNew Window

Fig. 3. Results of SMD simulation of GB1 and gpW proteins. (a) The force–extension curves of the unfolding process of protein GB1 show an force peak at the extension of ∼ 3.19 nm and forces all above 600 pN (710 pN, 714 pN, and 600 pN, respectively), which indicates the unfolding transition state. Three force–extension curves of the unfolding process of GB1 show very small difference in its unfolding process. (b) Six typical force–extension curves of the unfolding process of the protein gpW show multiple random force peaks. The peak forces are all below 200 pN, much smaller than that of GB1, which demonstrates that gpW is less mechanically resistant. In all SMD simulations, spring constant k = 100 kJ/mol/nm2 and pulling velocity v = 0.001 nm/ps.

Figure 3(b) shows several typical SMD simulation processes. In all simulations, the structure of protein gpW extends readily under the maximum stretching force of ∼ 200 pN, which is much smaller than the maximum unfolding force of protein GB1 (above ∼ 600 pN). Though in all simulations, a force smaller than 200 pN is enough to destroy the native structure of gpW gradually, but the detailed force–extension curves in each simulation show little similarities, which indicates that the unfolding processes of protein gpW have large variability. In contrast to gpW, the force peaks of GB1 with the same simulation parameters are always located at the extension of 3.19 ± 0.02 nm with maximal force of 660 ± 60 pN.%new simulations for GB1

To obtain a more statistic and quantitative understanding of the unfolding process of GB1, the free energy profile of GB1 unfolding was calculated through a series of umbrella sampling simulations. Neighboring simulations have enough overlapping to build PMF over extension range of ∼ 2.4 nm to ∼ 4.1 nm for GB1 and ∼ 1.3 nm to ∼ 4.5 nm for gpW.

The PMF profiles of GB1 and gpW are shown in Fig. 4. For GB1, PMF is very sensitive to extension. When the extension is about 4.1 nm, PMF is about 34 k_BT. On the other hand, PMF of gpW increases with force much more slowly. PMF of gpW is only 18 k_B T at the extension of 4.5 nm.

Fig. 4.

Figure Option ViewDownloadNew Window

Fig. 4. The PMF profiles of proteins GB1 and gpW along the coordinate of extension obtained from umbrella sampling. Typical histograms of conformations along the reaction coordinate of GB1 and gpW are shown in panels (a) and (b). Just a few umbrella sampling windows are selected to show a sufficient sampling along the coordinate. Totally 95 umbrella sampling simulations were performed for GB1 and 60 for gpW. (c) PMF of GB1 protein at zero force (red curve) shows a kink at extension of ∼ 3.18 nm, where the slope changes from ∼ 292 pN to ∼ 74 pN. At stretching forces of 30 pN, 60 pN, and 90 pN, PMF calculated from Eq. (4) is shown as blue, yellow, and gray curves, respectively. (d) PMF of gpW at zero force (red curve) shows a smooth curve which increases approximately linearly with extension. The slope of the curve is ∼ 29 pN. Blue, yellow, and gray curves show PMF at stretching forces of 10 pN, 20 pN, and 30 pN, respectively. Insets of both figures show typical conformations of GB1 and gpW with certain extension from the simulation.

The PMF of gpW is almost a linear function of extension when the extension is longer than 2 nm. The slope changes smoothly along extension from 1.4 nm to 4.5 nm. This smooth PMF of gpW is consistent with its downhill folding dynamics obtained in biochemistry study.

Different from gpW, the PMF profile of GB1 has a transition point at 3.2 nm where the slope changes from 292 pN to 74 pN, i.e., the slope drops to about one quarter of its value at extension smaller than 3.2 nm. This transition point is very close to the force peak position in the SMD simulations.

When a constant stretching force is applied to the end of protein in single molecule manipulation experiment or MD simulations, the free energy profile will tilt. Therefore, force-dependent free energy profile G_f(x) can be obtained from PMF (free energy profile at zero force, G₀(x)) according to the equation

The G_f(x) at different forces is also shown in Fig. 4. For GB1, the slope transition point becomes a peak when the force is 90 pN, which corresponds to the transition state of GB1 in its two-state folding/unfolding process.

As the PMF of gpW is almost a linear function of extension (Fig. 4(d)), its free energy profiles with stretching forces do not show any apparent barriers, which agrees with the downhill folding property of gpW. When the force is 20 pN, the PMF of gpW becomes almost a horizontal line. When the force is 30 pN or greater, the PMF decreases monotonically with extension. Three typical conformations with different extensions are shown for both proteins in the inset of Fig. 4.